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We show that arbitrary phase space vector fields can be used to generate 
phase functions whose ensemble averages give the thermodynamic tempera- 
ture. We describe conditions for the validity of these functions in periodic 
boundary systems and the Molecular Dynamics (MD) ensemble, and test 
G them with a short-ranged potential MD simulation. 
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I. INTRODUCTION 

The temperature of an equilibrium system is calculated from the mean kinetic energy of 
its particles. However, Rugh [] has recently derived a new expression whose average yields the 
temperature in the microcanonical ensemble. The derivation involves the differentiation of 
the phase space volume of the ensemble (whose logarithm gives the entropy) with respect to 
the energy. The resultant expression depends not only on the momenta of the particles, but 
also on their (spatial) configuration. This result reflects the fact that the temperature of a 
system effects the configurations adopted by that system, and consequently one would expect 
to be able to calculate the temperature from configurational as well as kinetic information. 

In [], a completely configurational version of Rugh's result was derived, and applied to 
Monte Carlo simulations, where the equipartition theorem cannot be used since only the 
configurational degrees of freedom are considered []. Not only did the temperature indeed 
correspond to the input temperature, but it was found to be a useful diagnostic for revealing 
coding errors as well. In [] , Rugh's result was applied to the non-equilibrium domain, as an 
alternative means of defining the local temperature. As such, it correctly accounted for heat 
fluxes where the kinetic temperature failed to do so. 

These applications have required extensions of Rugh's original work. In this paper we 
will justify these modifications. We will prove that the temperature expression used in [ 
is indeed equal to Rugh's expression to 0(1/N), and that the configurational temperature 
expression given by Eq. (8) in [] holds. 

For practical reasons, atomistic simulations are usually conducted under periodic bound- 
ary conditions. Furthermore, linear momentum is conserved in molecular dynamics (MD) 
simulations. Thus the ensembles explored by MD simulations are not the full microcanonical 
or canonical ensembles, but subsets of these (the "MD ensembles"). We will therefore also 
prove that the temperature expressions in [] hold for periodic boundary conditions, and in 
the canonical and microcanonical MD ensembles. 

In order to prove these results, we will derive a more general temperature expression 
in Sec. [H| Within the framework of this broader result, we will consider which functions 
yield the temperature for periodic systems and MD simulations. In Sec. |TT| we will use this 
theorem to obtain more specific results such as the equipartition theorem and Rugh's result. 
Finally, in Sec. |TV] we will apply these temperature- yielding functions to MD simulations of 
systems with short-range potentials. 

II. GENERALISED TEMPERATURE EXPRESSIONS 

Let us consider an iV-particle system at equilibrium. We define T = (Tx, . . . ,r 6 jv) = 
(pi, . . . ,P3N, Qii- ■ ■ i Q3n), where the g« and pi represent the 6N spatial coordinates and con- 
jugate momenta which determine the dynamics of the system via Hamilton's equations. 
The energy of our system is given by the Hamiltonian Ti.(T) = ^Zp\jm + V({qj}), where V 
represents the potential energy of the system. 

In this section we will prove the following result. Suppose we choose a vector field 
B(T) such that < | (VH • B(T)) | < oo, < | (V • B(T)) | < oo (where (. . .) represents 



an ensemble average) and (VT~C ■ B(T)) grows more slowly than e N in the thermodynamic 
limit. Then 

<™ • B ( r )) _ kT (1 ) 

(v-B(r)) ~ kT - (1) 

This result, as applied to the canonical ensemble, appears in [] without proof. For certain 
choices of B(T), it is possible to derive such temperature expressions using the approach 
of Gray & Gubbins [] - however, it is not evident how to extend their method of deriving 
hypervirial relations to arbitrary B(T). 

In Sec. [II A| and [II Q we will prove Eq. (|l[) for the microcanonical and canonical ensembles 



respectively. Moreover, it is possible to apply Eq. ([[]) to systems with periodic boundary 
conditions and the "MD ensembles" . We will develop the additional necessary conditions in 
Sec. |TTTJ and [TJJ. 



A. Microcanonical ensemble 

We begin with a proof of Eq. ([!]) in the microcanonical ensemble. Consider our A^-particle 
system, whose physical size is determined by a set of barriers or walls. If we denote by Q the 
set of all allowed T within our phase space, then the extent of Q in the spatial coordinates 
is limited by the physical size of the system. The momenta are unbounded, so that Q forms 
a cylinder in phase space. 

We define the surface of constant energy A(E) := {T : H(T) = E}, and the set of all 
points of equal or lower energy Q(E) := {T : Ti.(T) < E}. Thus Q = fi(oo). Traditionally, 
the microcanonical ensemble is the set of phase points whose energy lies between E and 
E + AE, AE <C E, with each point being equally likely to occur. Thus our (dimensionless) 
partition function would be 

■ dr. 



h 3N N\ J 

a(E+AE)\n(E) 

where h is Planck's constant. In the limit as AE — > 0, the microcanonical ensemble of 
energy E becomes the surface A(E) with a probability distribution given by the partition 
function 

1 r _. ._ 1 /■ cL4 



/ W( r)- B )dr = ^ / ^^_ / dw 



h 3N N ] J "V'-V / I h 3N N ] J ||VW(T) 

Q A(E) A(E) 

where dAg represents the infinitesimal area element on A(E) at F, and d/i£ = 
cL4g;/||V?i(r)||. We note immediately the dimensional inconsistency of this expression. 
From a physical point of view, it is necessary at this point to reduce the spatial coordi- 
nates and momenta to dimensionless units p* and q*. Consequently, all functions of these 
coordinates will remain dimensionless, and can be converted to the appropriate units after 
calculation in reduced units. Physically, this is equivalent to selecting a standard basis of 
units by which to measure all the Pi, and all the qi (and consequently, the units of our 



resultant temperature). However, this choice of units for the pi and qi can always be made 
independently of each other, so that this does not present a problem in practice. Indeed, 
since there is no unique choice of units for converting the phase space coordinates to a 
reduced form, there will be infinitely many different instantaneous expressions for the tem- 
perature which will all give equivalent values in different bases of units (see, e.g., []). In 
what follows, we will assume that our pi and qi are dimensionless as required (as will be the 
factor h 3N ). 

In the surface ensemble, the entropy S(E) can be defined as 



,3{E)/k 



h 3N N\ 



dfi E , 



A(E) 



and the average value of a phase function B(T) in the ensemble is 

J A(E) B(T)du E e -s(*)A 



(BiT)), 



Ia(E) d^E 



h m N\ 



B(T)dfi E . 



A(E) 



The temperature of a microcanonical ensemble with energy E, temperature T, entropy S, 
and volume V is determined via the thermodynamic relation 



dS(E) 



T(E) dE 

Suppose, for an arbitrary vector field B(r), we define B(T) = V7Y(r) • B(r) and 

1 



W B (E) 



,S B {E)/k 



h 3N N\ 



B(T) dfi E , 



(2) 



(3) 



A(E) 



where we assume < (B(T)) E < oo. From first principles, 



h m N\ 



dW B {E) 
dE 



lim h m N\ 

5^0 



lim - 

5^0 S 



W B (E + 5)- W B {E) 



lim- 

5^0 S 



B(T) dfi E+5 - J B(T) dfi E 

A(E+S) A(E) 



B(r) • n(r) dA E+s - J B(r) • n(r) dA E 

A(E+6) A(E) 



where n(r) is a unit normal vector to the surface A(E) at I\ We apply Gauss' Theorem to 
obtain 

h™N\ ^^ = lim i T f V • B(T) dA, dC 
dE s^o 6 J J ? 

E A(0 

= J v-B(r)d/i £ . 

A(E) 



Therefore it follows that 

1 _ ldS B (E) _ (V-B(r)) J 



kT B (E) ' k dE (VH ■ B(r))j 



(4) 



From Eq. @, we have that S B (£) = S(E) + kin (B(T)) E . Now, as long as {B(T)) E grows 
more slowly than e N in the thermodynamic limit, dS B (E)/dE = dS(E)/dE in the thermo- 
dynamic limit (a relation we will henceforth denote as dS B (E)/dE ~ dS(E)/dE). Conse- 
quently, Tjg(E) ~ T(E), and we recover the same result as Eq. ([Xj) . Furthermore, we may 
drop the condition that (B(T)) E be positive, since Eq. (f|) gives the same value, whether we 
use ±B(r). 

In order to apply Gauss' Theorem, we require that VTL be continuous (ie that TC be 
differentiate) for finite energies. As far as conditions on B(T) are concerned, we require that 
| (V • B(T)) E | < oo, that < | (VH ■ B(T)) E | < oo, and that (In | (\7H ■ B(T)) \)/N ~ 0. 
Given the dependence of A(E) on the Hamiltonian, the family of vector fields B(T) that 
obey these criteria is not immediately obvious, and we have not attempted to develop a 
more general method of generating such vector fields. We simply note that the last condition 
allows all finite-order polynomials and bounded functions of the Pi and g,, as well as ratios 
of finite-order polynomials. However, it remains clearer from the canonical case (where the 
domain of integration does not depend on the Hamiltonian) whether such functions will 
obey the first two conditions. 

B. Microcanonical Periodic and MD Systems 

Let us now consider the necessary changes to the proof of Eq. (|j) in order for it to hold 
in a periodic system. For periodic systems, the extent of Q in the spatial coordinates is no 
longer determined by boundary walls, but by the size and shape of the primitive cell. If the 
primitive cell of the periodic system is the same size and shape as the bounded system, then 
fl will be the same in both cases. 

The difference between the bounded and periodic systems is that particles cannot pass 
through the walls of the bounded system. This implies that the energy at the walls is infinite, 
so that our surfaces of constant energy lie entirely within fl, and do not pass through the 
boundary, which we denote dfl. This assumption is implicit in our application of Gauss' 
Theorem. 

In the periodic system, particles can pass through the "walls" of our primitive cell, 
reappearing on the other side of the cell. Therefore surfaces of constant energy can (and 
do) pass through dfl. Thus, when we use Gauss' Theorem, our Gaussian surface consists 
not only of A(E) and A(E + h), but also all points on dQ whose energies lie between E and 
E + h. This extra term is of the form 



lim — 
h-*o h 



f B(r) • n(r) cL4 - f B(r)-n(r)cL4 = f B(r)-n(r)dL 



dU(E+h) dQ(E) 



dA{E) 



where dQ(E) = d£l fl £l(E), dA is the volume measure on <9f2, dA(E) = d£l fl A(E), and 
dL is the volume measure on dA(E). Note that n(r) does not necessarily point in the same 



direction as V7i(T), since the walls of the primitive cell are not determined by the energy 
surfaces. For Eq. (|4]) to hold in every microcanonical ensemble, we require, as a condition 
on B(T), that 



B(r) ■ n(r) AL = ME. (5) 

8A(E) 

To determine which functions satisfy this criterion, we consider a system where one of the 
particles is at one of the walls of the primitive cell [], corresponding to a phase point T a . 
There is an equivalent system where this particle is placed on the "opposite" wall of the 
primitive cell, represented by F b . It follows that n(T a ) = — n(T 6 ). Therefore, since T a and 
r^ must lie in the same microcanonical ensemble, if B(T a ) = B(TJ, then the criterion of 
Eq. (^) is satisfied. Therefore any function which is periodic in the primitive cell [ie such 
that, if T a and T b describe the same state, then B(T a ) = B(T fe )] will satisfy Eq. ([[]) for 
periodic systems. Note that this is a sufficient condition but not a necessary one. 

Finally, let us consider the MD microcanonical ensemble. This ensemble represents 
the family of systems encountered during a constant energy molecular dynamics simulation, 
where linear momentum is conserved. Let Omd be the set of allowed T for such a simulation. 
Clearly Qud is smaller than Q, which admits all possibilities for the total linear momentum. 
Phase points T with the same linear momentum in the ^-direction, say, all lie on the same 
(hyper)plane in Q, so that Qmd is the intersection of Q with the three phase space planes 
which correspond to conservation of linear momentum in each Cartesian direction. We 
denote their normal vectors as P x , P y , and P z . 

The entropy will still correspond to the phase space volume, except that this volume is 
now QN — 4 dimensional. However, we can only apply Gauss' Theorem to the projection of 
the vector field B(T) onto I^md- Alternatively we must select B(T) so that it lies entirely 
in Omd- Such a B(T) must satisfy the condition that B(T) • P a = (a = x, y, z). In this 
case, Eq. (|J) will generate the correct temperature in the MD ensemble. 

C. Canonical Ensemble 

We now move on to a proof of Eq. (|l|) in the canonical ensemble, starting with the 
bounded case. We invoke Gauss' Theorem over Q(E) for an arbitrary vector field in phase 
space B(r)e~ /3W( ^ r ^ (where we assume a finite, positive /?), ie 



J e -^( r )B(r) • n(r) 6A E = J v ■ (B(r)e~^( r )) dr 

A(E) Cl(E) 

= J e -^ r V • B(r) dr - (3 J e-^^vrUT) ■ B(r) dr. 

Q(E) O(E) 

In the limit as E — » oo, we obtain 
Jirn^ e- pE J B(T) • n(r) dA E = J e' m{T) V ■ B(T) dT - (3 f e -^ r VW(r) • B(T) dr. 

A(E) Q Q 

(6) 



For Eq. (|6|) to be of any use, we require that the two integrals on the right hand side be 
finite. For the latter integral this gives us 

/ e -/*w(r)vft(r).B(r)dr <oo =► lim / e~ 0E B(T) • n(r) 6A E = 0. 

J E^co J 

n a(_b) 

This means that whenever the last integral in Eq. @ exists, the left hand side of Eq. (||) 
must be identically zero. It follows by rearrangement that 

JL B n /^'v-B(r)dr (v . B(r)> 

kr p /e-^( r )v^(r)-B(r)dr " (vw(r) • B(r)) ' lJ 

in agreement with Eq. (|]). Since /3 is finite, we have subsumed the first two conditions on 
B(r) into the proof. The third is implicit in the convergence of (VTC ■ "B(T)) E . For Eq. (|7|) 
to hold, we require that the integral 

e -PE J B(r) . fi(]?) AAe = e P[TS { E)-E] {B{t))e = e P[TS B (E)-E] 
A(E) 

converge. When we consider that e^ TS ^~ E ^ also converges, but that e s ^ does not, we 
immediately obtain the third condition on B(T). Thus the conditions for Eq. ([]]) to hold in 
the canonical ensemble are the same as those for the microcanonical ensemble. 



D. Canonical Periodic and MD Ensembles 

As with the microcanonical case, we must be careful in our application of Gauss' Theorem 
to canonical systems with periodic boundary conditions. In analogy with the microcanonical 
case, our Gaussian surface consists not only of A(E), but of dQ(E) as well, and the left hand 
side of Eq. (|j) becomes 

lim / e~ m ^B(T) ■ h(T) dA E . 

A(E)uan(E) 

We have already seen that the integral over A(E) must go to zero in order for (V7i(T) ■ B(T)) 
to exist, so we simply require that 

J e -^( r )B(r) ■ n(r) dA = V/3. 
an 

However, via the properties of the Laplace transform we have that 

J e" /3W(r) B(r) ■ n(r) dA = V/3 <£> I B(T) • n(T) dL = ME. 

dQ dA(E) 



Therefore the condition under which Eq. (P will hold in all canonical ensembles is equivalent 
to the condition under which Eq. ([!]) will hold in all microcanonical ensembles. Just as in the 
microcanonical case, Eq. (|7|) will hold in the canonical ensemble as long as B(T) is periodic 
in Q. 

Finally, we consider the canonical MD ensemble. As with the microcanonical MD ensem- 
ble, our application of Gauss' Theorem requires that B(T) lie in Qmd, so we again require 
that B(r) • P Q = (a = x, y, z) for Eq. ([?p to hold in the canonical MD ensemble. 

Thus the conditions for Eq. ([[]) to hold for the periodic boundary system and the "MD 
ensembles" are the same in both the canonical and microcanonical ensembles. 

III. FORMULAE 

Having proven Eq. (p]) in the canonical and microcanonical ensembles, and found the 
conditions for it to hold in systems with periodic boundary conditions and the MD ensembles, 
we now demonstrate its use in generating expressions whose phase space average yields the 
system temperature. 

If we choose B(T) = (0, . . . , Tj, . . . , 0), so that only the i-th component is non-zero, then 
we obtain 

_ (vw(r)-B(r)) _ / m\ 
equip " (v-B(r)) 'y'dvj [) 

This is the familiar Generalised Equipartition Theorem. If Tj is a momentum, then we obtain 
the Equipartition Theorem, (pj/m) = kT. If it is a coordinate, then we obtain the lesser 
known Clausius Virial Theorem, (— ^Fj) = kT, where Fi is the generalised force acting on 
coordinate g« []. We note that the Clausius Virial Theorem gives a function of coordinates 
only, whose average is the temperature of the system. However, the function B(T) (= T) is 
not periodic in Q in this case, so that this theorem does not hold for periodic systems. It is 
therefore of little use to practitioners of most MD simulations as a means of calculating the 
temperature. 

If we select an arbitrary vector field X(T), and choose 

B(r) = vw(r)-x(r)' (9) 

then for all choices of X(T), B{T) = 1. Consequently, we obtain 

1 / X(T) \ 

vf = V vw(r)-x(r)/' (10) 

providing this average exists. Substituting X(T) = V7Y(T), we obtain Rugh's final equation 
[]. Since the Hamiltonian is periodic in systems with periodic boundary conditions, B(T) will 
also be periodic, so that Rugh's result holds in periodic systems. Furthermore, it satisfies 
the criterion for the MD ensembles, so that it can be applied to MD simulations as well. 



IV. EXAMPLE : SIMULATION APPLICATION 

In this section we consider the application of Eq. (|I|) to a simulation of a system of 
particles interacting with a short-range pair potential, as in []. These simulations employ 
periodic boundary conditions, and as a consequence, the forces acting on a body are not 
correlated with the absolute positions of the particles, but only their relative positions. Thus 
many of the simple vector fields whose divergences are easily calculated (such as B(T) = T) 
do not satisfy the first criterion for Eq. (|1|) to hold; in this case, (VTC ■ T) E = 0. In general, it 
is more difficult to find functions which are correlated with the inter-particle forces. Due to 
this difficulty, from this point on we restrict ourselves to choices of B(T) which are directly 
related to VH, to ensure that this condition is met. 

A. Theory 

In systems of particles interacting with a short-range pair potential $(r), the Hamilto- 
nian can be separated into a momentum contribution (the kinetic energy K) and a spatial 
contribution (the potential energy V), ie 

37V 2 TV 

W(T) = K({ Pi }) + V({q t }) = £ g- + ££ *(M). (11) 

where r^ = r* — Vj, and r, is the vector describing the position of the i-th particle. If 
the potential has a continuous first derivative, then VH satisfies the requirements of Gauss' 
Theorem, and consequently those of our temperature expressions. Note that if we define r y 
as the minimum image separation of the i-th and j-th particles, then V is periodic in the 
spatial coordinates. Thus VTL will be periodic as well. As a consequence, we can obtain the 
temperature of our system using Rugh's expression, ie by substituting X(r) = VH(T) into 
Eq. (|TDD above. 

We now make the following important observation — since V7i satisfies the criteria for 
Eq. ([!]) to hold in periodic boundary systems and MD ensembles, it follows that VK and 
VV must as well. Therefore, we would expect to be able to generate the temperature by 
substituting X(r) = VK({ Pi }) and X(r) = W({ft}) into Eq. (0). 

Since the interaction potential is short-ranged, V will grow as iV in the thermodynamic 
limit. Consequently, the Hamiltonian grows as N in the thermodynamic limit, and if we 
substitute B(r) = VH(T), B(T) = VK({ Pi }) and B(r) = W({&}) into Eq. (g), we would 
also expect to generate the temperature. 

In this paper we will not examine the temperatures generated from the kinetic energy, 
since they are closely related to the equipartition temperature [Eq. (§)], and do not reveal 
any new results. Our interest lies in the fact that temperature expressions generated with 
W({qi}) contain no explicit reference to the momenta in our system, a fact which has 
been exploited in []. The temperature we obtain from substituting X(r) = V7Y(r) into Eq. 
( [TOD we denote by T nor R - - "nor" since it is generated using the normal vector field VTt, 
and "R" since it is generated using Rugh's prescription. In a similar manner, we denote by 
^conR the temperature we obtain from substituting X(r) = W({gj}) — the configurational 



part of the Hamiltonian -- into Eq. fllPl). When substituting these vector fields into Eq. 
(H), we denote the corresponding temperatures as T norF and T conF , "F" denoting that we are 
calculating a ratio (fraction) of averages in this case. 

In making the appropriate substitutions, we obtain the following expressions : 

(12a) 



(12b) 
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/-E,V 4 - 
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KJ- conR 
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1 _( 


f-E.V.-F,) 








n± norF 


\ E« ^T + Fj ) 





(12c) 



EiV4 ' Fl) , (12d) 



fcT conF (£,F 2 ) ' 

where the label % refers to the particle, rather than the generalised coordinate. Fj represents 
the (vector) force acting on particle i, p, represents its momentum, V» = [-£r, ^-, ^-], where 
Xi,yi and Zi refer to the Cartesian coordinates of r i; and : represents the dyadic operator 
(ie for vectors a, b and matrix M, ab : M = Ea a a a^B^-0a) ■ Eq. ( |12d| ) corresponds to the 
temperature expression used in [], and Eq. (|12cj) corresponds to the temperature expression 
used in []. 

If we consider the second term on the right hand side of Eqs. ( |12ayi2b|) , the numerator 
increases as N for a short-ranged potential (since FjFj : VjFj will not contribute anything 
at large particle separations), but the denominator increases as N 2 . Therefore this second 
term becomes negligible in the thermodynamic limit. Thus the order 1 term is contained in 
the first term on the right hand side of Eqs. ( p.2a|J12b| ). We will denote by T norl and T conl 



the temperature calculated by the omission of these second terms respectively, ie 

1 /^-E.V.-FA 



kT nml \ E,F 2 + ^ 
1 /-E.V.-F 



*■*■ conl \ 2—ii " i 



(14e) 



(14f) 



We expect that the temperatures given by Eqs. (|12a|) -( |l4~i]) should all be equal in the 



thermodynamic limit. It is therefore of interest to compare their rates of convergence to this 
limit, in order to ascertain the appropriateness of their use. 



10 



B. Results 

As an application of the above theory, we considered a three-dimensional microcanonical 
WCA-potential system. The WCA potential is defined as follows [], 



4e 


'/ \ 12 


- (f)' 


+ e, r < 2 1 /V 


o, 






otherwise 



$(r) = I " C [W W J TC ' , (13) 

[ 0, otherwise 

where a and e represent our units of length and energy respectively. This potential is 
continuous, has a continuous first derivative and a piecewise continuous second derivative. 
Due to the discontinuity in the derivative of the force, errors appear in the computed system 
trajectories whenever the separation between two particles crosses the r = 2 1//6 <r boundary. 
However, these errors are too small, in comparison with system size errors, to affect our 
results. Thus, if we substitute this pair potential into Eq. (|TT|), then we expect each of the 
temperatures defined in Eqs. ( |12a| )- (|14i|) to be equal, to order (In N)/N. 



Values of these six temperatures were calculated for the 3D microcanonical simulation of 
a periodic WCA system at various sizes N, reduced number densities p*, and reduced total 
energies per particle E*. They were determined by the average of ten separate simulations, 
each of 200000 timesteps (of 8t* = 0.001). The errors associated with each temperature were 
given by one third of the maximum deviation from the average over these ten runs. 

The first comparison was made between systems with the same density and size, but 
differing energies. The values of the equipartition and normal temperatures (T equip , T nor R, 
T nor F and T nor i) were calculated for a system of 500 particles with reduced density p* = 0.8, 
and reduced energies per particle ranging from E* = 0.8 to E* = 2.5. These values appear 
in Table [II The four temperatures agree to within 0.6-0.8% of the equipartition temperature 
over the range of energies shown. 

The values for the three configurational temperatures match the corresponding normal 
temperatures to within 0.01%, ie to the number of digits shown in Table Q. This can be 
explained in terms of the kinetic and configuration terms in the numerator and denominator 
of the normal temperature expressions. At high densities, the configuration terms are much 
larger than those contributed by the momentum terms — in two dimensions, this is typically 
a difference of four orders of magnitude, and in three dimensions, the difference is about 
six orders of magnitude. For this reason, the value of the normal temperature can be 
considered as a "perturbation" to the corresponding configurational temperature which has 
a negligible effect on our results. It is interesting to note, given this dependence on the 
physical structure of the system rather than on its momentum distribution, that the normal 
and configurational temperature expressions yield the correct temperature across the solid- 
liquid phase transition, despite the difference in the microscopic arrangements of atoms on 
either side of the transition temperature. 

In Fig. [3], we compare a series of systems of fixed energy per particle (E* = 1.5) and 
system size (N = 864), but with varying densities. The discrepancy between T norl and 
T nor R increases when the density of the system is decreased — while T nor R and T nor F agree 
quite well with the equipartition values, T nor i becomes less and less reliable. However, 
in the thermodynamic limit, T nor i must converge to the other two normal temperatures. 
This result indicates that, while T norl and T norR must converge towards the thermodynamic 

11 



temperature, irrespective of the density, larger systems sizes are required for the same degree 
of convergence of T norl as the density drops. 

We should also note from Fig. [J] that, while T norR and T norF are indistinguishable from 
their configurational counterparts on the scale of the graph (and hence are not shown), the 
difference between T nor i and T con i becomes evident below densities of p* ~ 0.5. This is a 
result of the drop in the number of particle interactions per timestep at lower densities. 
When the number of these interactions is reduced, the configurational contributions do not 
dominate the kinetic contributions as they do in the high density regime. Consequently, the 
inclusion of kinetic terms (which, by themselves would produce a value within 0.1% of the 
equipartition value) in T norl will always correct T conl towards the equipartition value. 

To further examine the system size dependence of our temperature expressions, we con- 
sider a single state point (p* = 0.8, E* = 1.5), and compare the temperature expressions as a 
function of the number of particles in the system, ranging from N = 108 to iV = 2048. The 
results of this comparison appear in Fig. |^, where the three configurational temperatures 
have been plotted against inverse system size. At this density, the difference between the 
normal temperatures and the corresponding configurational temperatures is not distinguish- 
able on the scale of the graph for all but the 108 particle system (where the discrepancy is 
0.02%), so we show only the configurational temperatures. We observe, within the errors 
of our calculations, the convergence of all four temperatures towards a common value. We 
would interpret this value as the thermodynamic temperature of a system at that state 
point, in the thermodynamic limit. 

V. CONCLUSION 

We have derived a general functional which, given a vector field B(T) which satisfies 
certain broad conditions, will determine the thermodynamic temperature of an equilibrium 
system in the thermodynamic limit via Eq. (H). Its rate of convergence in the thermodynamic 
limit will be determined by the order of (B(T) • VH). We note, however, that if we define 
B(T) as per Eq. (||), then (B(T) • V7Y) = 1, and what we obtain in Eq. (P is precisely the 
derivative of the logarithm of the ensemble phase space volume with respect to the energy. In 
the thermodynamic limit, this will yield the thermodynamic temperature dS/dE. However, 
for different B(T), the value we obtain will depend upon our sampling of phase space during 
the simulation, and hence the values obtained from different expressions may vary. The 
temperature expressions T norR and T conR fall into this category. 

One practical problem that arises from the application of Eq. ([J) to periodic boundary 
systems is the difficulty in avoiding vector fields B(T) such that (B(T) • VH) = 0. To 
circumvent this problem, we have only considered vector fields B(T) that are linear trans- 
formations of VH. This approach is by no means exhaustive, but serves to demonstrate one 
application of this theory. 

It is clear from Eqs. ( |12a| )- (|14fl) that T norR will be computationally more expensive than 



Tnori or TnorF - _ the omitted term involves calculations which assume the intermolecular 
forces to have already been evaluated, thus requiring a second force loop. It is therefore of 
interest to determine whether these approximations to dS/dE make a useful substitution for 
T nor R. From our results we conclude that T norF is more reliable than T norl , and in our work 



12 



is a useful expression for the temperature whenever T nor R is valid. It is for these reasons 
that the fractional forms (T norF or T conF ) appear in []. 

Eq. ([!]) has important consequences for practitioners of non-equilibrium MD simulations, 
stressing the fact that the instantaneous kinetic energy per kinetic degree of freedom is not 
the only function whose ensemble average yields the temperature. The preference given 
to the kinetic energy is generally due to its ease of calculation: apart from this, there is 
no reason --in both equilibrium and non-equilibrium calculations --to prefer the kinetic 
energy expression over any other. 
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TABLES 

TABLE I. A comparison of values of the three normal temperatures with values of the equipar- 
tition temperature, for simulations of systems with size N = 500, reduced density p* = 0.8, and 
various reduced energies per particle E* . For each normal temperature, two values are reported. 
The first is the temperature as determined from the simulation, and the second is the discrepancy 
between that normal temperature and the equipartition temperature, given as a percentage of the 
equipartition temperature. The numbers in brackets indicate the error in the last decimal place 
given. 
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FIG. 1. Variation of temperature values with system density. For a system of 864 particles at 
reduced energy per particle E* = 1.5, various temperatures are given for different reduced densities. 
Temperatures are reported as a fraction of the equipartition temperature. 
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FIG. 2. Variation of temperature values with system size. For the state point p* = 0.8, 
E* = 1.5, the three configurational temperatures and the equipartition temperature are given 
for different system sizes. 
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